Studies 1-4: Figures

Luhrmann, Weisman, et al.

In this notebook, we create figures for Studies 1-4.

source("./scripts_general/dependencies.R")
source("./scripts_general/custom_funs.R")
source("./scripts_general/var_recode_contrast.R")
source("./study1/scripts_s1/s1_var_groups.R")
source("./study2/scripts_s2/s2_var_groups.R")
source("./study3/scripts_s3/s3_var_groups.R")
source("./study4/scripts_s4/s4_var_groups.R")

setwd("./study1/analysis/")

The working directory was changed to /Users/anthrouser/Documents/Projects/Mind and Spirit overview paper/mindspiritquant/study1/analysis inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

source("../../scripts_general/data_load.R")

Parsed with column specification:
cols(
  .default = col_double(),
  date = col_date(format = ""),
  researcher = col_character(),
  country = col_character(),
  site = col_character(),
  religion = col_character(),
  subject_gender = col_character(),
  subject_job = col_character(),
  subject_schedule = col_character(),
  subject_livedhere = col_character(),
  subject_lang = col_character(),
  subject_marital = col_character(),
  subject_hs = col_character(),
  subject_liveswith = col_character(),
  servicesperweek = col_character(),
  study = col_character()
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
  .default = col_double(),
  date = col_date(format = ""),
  researcher = col_character(),
  country = col_character(),
  quad = col_character(),
  subject_gender = col_character(),
  subject_job = col_character(),
  subject_schedule = col_character(),
  subject_livedhere = col_character(),
  subject_lang = col_character(),
  subject_marital = col_character(),
  subject_hs = col_character(),
  subject_school = col_character(),
  subject_liveswith = col_character(),
  servicesperweek = col_character(),
  godexpviaawe_freq = col_character(),
  sleephabit = col_character(),
  researcher_date = col_date(format = ""),
  researcher_notes = col_character(),
  site = col_character(),
  religion = col_character()
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
  .default = col_double(),
  ctry = col_character(),
  wher = col_character(),
  recr = col_character(),
  whoc = col_character(),
  demo_chur = col_character(),
  demo_ethn = col_character(),
  demo_maj = col_character(),
  demo_pocc = col_character(),
  demo_rlgn = col_character(),
  demo_sex = col_character()
)
See spec(...) for full column specifications.
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  .default = col_character(),
  X1 = col_double(),
  epi_subj = col_double(),
  epi_demo_age = col_double(),
  epi_demo_ses_num = col_double(),
  epi_demo_howr_num = col_double(),
  epi_demo_affr_num = col_double(),
  epi_demo_tung_num = col_double(),
  epi_demo_ytng_num = col_double(),
  epi_demo_affr_cat = col_logical(),
  epi_demo_tung_cat = col_logical(),
  epi_demo_ytng_cat = col_logical(),
  epi_version = col_double(),
  epi_charc = col_double()
)
See spec(...) for full column specifications.
2 parsing failures.
 row       col expected  actual                           file
1025 epi_charc a double Unclear '../../study3/data/d_demo.csv'
1027 epi_charc a double Unclear '../../study3/data/d_demo.csv'
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  epi_ctry = col_character(),
  epi_subj = col_double(),
  score = col_double()
)
Joining, by = c("epi_ctry", "epi_subj")
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  epi_ctry = col_character(),
  epi_subj = col_double(),
  score = col_double()
)
Joining, by = c("epi_ctry", "epi_subj")
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  epi_ctry = col_character(),
  epi_subj = col_double(),
  question = col_character(),
  response = col_double(),
  order = col_double(),
  question_text = col_character()
)
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  epi_ctry = col_character(),
  epi_subj = col_double(),
  question = col_character(),
  response = col_double(),
  order = col_double(),
  question_text = col_character()
)
Joining, by = c("epi_ctry", "epi_subj", "question", "response", "order", "question_text")
Joining, by = c("epi_ctry", "epi_subj")
Column `epi_ctry` joining character vector and factor, coercing into character vectorJoining, by = "epi_subj"
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  .default = col_double(),
  p7_ctry = col_character(),
  p7_abs_check = col_character(),
  p7_dse_check = col_character(),
  p7_se_check = col_character(),
  p7_unev_check = col_character(),
  p7_exsen_check = col_character(),
  p7_por_check = col_character(),
  p7_mm_check = col_character(),
  p7_dem_sex = col_character(),
  p7_dem_pocc = col_character(),
  p7_dem_major = col_character(),
  p7_dem_ethnicity = col_character(),
  p7_dem_rur.urb = col_character(),
  p7_dem_affrd.basics = col_character(),
  p7_dem_religion = col_character(),
  p7_dem_church = col_character(),
  p7_dem_holy.tung.gif = col_character(),
  p7_abs_child.exp_cat = col_logical(),
  p7_abs_poetic_cat = col_logical(),
  p7_abs_tv.real_cat = col_logical()
  # ... with 162 more columns
)
See spec(...) for full column specifications.
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("study", "p7_ctry", "p7_subj", "abs_score", "cog_score", "ctl_score", "dse_score", "hall_score", "para_score", "por_score", "pv_score", "spev_score")
Joining, by = c("study", "p7_ctry", "p7_subj", "abs_score", "cog_score", "ctl_score", "dse_score", "hall_score", "para_score", "por_score", "pv_score", "spev_score")

d_all <- d1 %>%
  select(study, country, site, religion, subject_id, 
         pv_score, abs_score, spev_score) %>%
  mutate(religion = recode_factor(religion,
                                  "indigenous" = "Indigenous Religion",
                                  "charismatic" = "Charismatic Christianity"),
         # rescale to 0-1
         pv_score = pv_score/3) %>%
  full_join(d2 %>% 
              select(study, country, subj, 
                     abs_score, spev_score, dse_score) %>% 
              rename(subject_id = subj) %>%
              # rescale to 0-1
              mutate(spev_score = spev_score/4,
                     dse_score = dse_score/5,
                     religion = "General Population")) %>%
  full_join(d3 %>% 
              select(study, epi_ctry, epi_sample, epi_subj, 
                     por_score, spirit_score) %>%
              rename(country = epi_ctry, 
                     religion = epi_sample,
                     subject_id = epi_subj,
                     spev_score = spirit_score) %>%
              mutate(religion = recode_factor(religion,
                                              "general population" = "General Population",
                                              "charismatic" = "Charismatic Christianity"))) %>%
  full_join(d4 %>%
              select(study, p7_ctry, p7_subj, 
                     por_score, pv_score, abs_score, spev_score, dse_score) %>%
              rename(country = p7_ctry,
                     subject_id = p7_subj) %>%
              # rescale to 0-1
              mutate(por_score = por_score/2,
                     pv_score = pv_score/3,
                     spev_score = spev_score/4,
                     dse_score = dse_score/5,
                     religion = "General Population")) %>%
  mutate(religion = factor(religion, 
                           levels = c("General Population", 
                                      "Indigenous Religion",
                                      "Charismatic Christianity")),
         study = gsub("study", "Study", study)) 

Joining, by = c("study", "country", "religion", "subject_id", "abs_score", "spev_score")
Column `religion` joining factor and character vector, coercing into character vectorJoining, by = c("study", "country", "religion", "subject_id", "spev_score")
Column `religion` joining character vector and factor, coercing into character vectorJoining, by = c("study", "country", "religion", "subject_id", "pv_score", "abs_score", "spev_score", "dse_score", "por_score")

Summary statistics

By study

Study 1

plot_s4_spev <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = country, y = spev_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(spev_score, na.rm = T),
                              sd = sd(spev_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 1), breaks = seq(0, 1, 0.25)) +
  labs(x = "Country", y = "Spiritual Events")

plot_s4_dse <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = country, y = dse_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(dse_score, na.rm = T),
                              sd = sd(dse_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 1), breaks = seq(0, 1, 0.25)) +
  labs(x = "Country", y = "Daily Spiritual Experiences")

plot_s4_pv <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = country, y = pv_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(pv_score, na.rm = T),
                              sd = sd(pv_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 1), breaks = seq(0, 1, 0.25)) +
  labs(x = "Country", y = "Porosity Vignettes")

plot_s4_por <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = country, y = por_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(por_score, na.rm = T),
                              sd = sd(por_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 1), breaks = seq(0, 1, 0.25)) +
  labs(x = "Country", y = "Porosity Scale")

Spiritual Events scores

plot_s4_abs <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = country, y = abs_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(abs_score, na.rm = T),
                              sd = sd(abs_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 1), breaks = seq(0, 1, 0.25)) +
  labs(x = "Country", y = "Absorption")

plot_grid(plot_s4_spev, plot_s4_dse, plot_s4_pv, plot_s4_por, plot_s4_abs, NULL, ncol = 2, labels = c("A", "B", "C", "D", "E"))

Removed 2 rows containing missing values (geom_point).

Study 1

d1 %>%
  ggplot(aes(x = country, y = spirit_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(spirit_score, na.rm = T),
                              sd = sd(spirit_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, 0.25)) +
  labs(x = "Country", y = "Study 1: Spiritual Events score (range: 0-1)",
       caption = "Error bars are ±1 standard deviation from the mean")

Study 2

d2 %>%
  ggplot(aes(x = country, y = spev_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(spev_score, na.rm = T),
                              sd = sd(spev_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 4), breaks = seq(0, 4, 1)) +
  labs(x = "Country", y = "Study 2: Spiritual Events score (range: 0-4)",
       caption = "Error bars are ±1 standard deviation from the mean")

Study 3

d3 %>%
  # correct for scaling in original dataset
  mutate(spirit_score = spirit_score * 4) %>%
  ggplot(aes(x = epi_ctry, y = spirit_score, color = epi_ctry, fill = epi_ctry)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(epi_ctry) %>%
                    summarise(mean = mean(spirit_score, na.rm = T),
                              sd = sd(spirit_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 4), breaks = seq(0, 4, 1)) +
  labs(x = "Country", y = "Study 3: Spiritual Events score (range: 0-4)",
       caption = "Error bars are ±1 standard deviation from the mean")

Study 4

d4 %>%
  ggplot(aes(x = p7_ctry, y = spev_score, color = p7_ctry, fill = p7_ctry)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(p7_ctry) %>%
                    summarise(mean = mean(spev_score, na.rm = T),
                              sd = sd(spev_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 4), breaks = seq(0, 4, 1)) +
  labs(x = "Country", y = "Study 4: Spiritual Events score (range: 0-4)",
       caption = "Error bars are ±1 standard deviation from the mean")

All studies

d_all %>%
  ggplot(aes(x = country, y = spev_score, color = country, fill = country,
             group = study)) +
  geom_point(position = position_jitterdodge(jitter.width = 0.25,
                                             jitter.height = 0,
                                             dodge.width = 0.75), 
             alpha = 0.1) +
  geom_pointrange(data = d_all %>%
                    group_by(study, country) %>%
                    summarise(mean = mean(spev_score, na.rm = T),
                              sd = sd(spev_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd, shape = study), 
                  position = position_dodge(width = 0.75),
                  color = "black") +
  scale_shape_manual(values = 21:24) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  guides(color = F, fill = F, 
         shape = guide_legend(override.aes = list(fill = "black"))) +
  theme(legend.position = "bottom") +
  labs(x = "Country", y = "Spiritual Events score (rescaled to 0-1)",
       shape = "Study", 
       caption = "Error bars are ±1 standard deviation from the mean")

d_spev_sum <- d_all %>%
  filter(!is.na(spev_score)) %>%
  group_by(study, country, religion) %>%
  summarise(mean = mean(spev_score, na.rm = T),
            sd = sd(spev_score, na.rm = T),
            n = n()) %>%
  ungroup()

d_all %>%
  ggplot(aes(x = country, y = spev_score, 
             color = country, fill = country,
             group = religion)) +
  facet_grid(~ study, scales = "free", space = "free") +
  geom_point(position = position_jitterdodge(jitter.width = 0.8,
                                             jitter.height = 0.02,
                                             dodge.width = 0.75), 
             alpha = 0.15) +
  geom_pointrange(data = d_spev_sum,
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd, 
                      shape = religion), 
                  position = position_dodge(width = 0.75),
                  fill = "black",
                  color = "black") +
  geom_text(data = d_spev_sum %>%
              mutate(ypos = case_when(
                grepl("charismatic", tolower(religion)) ~ mean + sd + 0.05,
                TRUE ~ mean - sd - 0.05)),
            aes(y = ypos, label = paste0("n=", n)), 
            position = position_dodge(width = 0.75),
            size = 3, color = "black") +
  scale_shape_manual(values = 21:24) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  guides(color = F, fill = F, 
         shape = guide_legend(override.aes = list(fill = "black"))) +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(x = "Country", y = "Spiritual Events score (rescaled to 0-1)",
       # caption = "Error bars are ±1 standard deviation from the mean",
       shape = "Religion")

Porosity scores

d_por_sum <- d_all %>%
  gather(por_scale, por_score, c(pv_score, por_score)) %>%
  mutate(por_scale = recode(por_scale,
                            "pv_score" = "Porosity Vignettes",
                            "por_score" = "Porosity Scale")) %>%
  filter(!is.na(por_score)) %>%
  group_by(study, country, religion, por_scale) %>%
  summarise(mean = mean(por_score, na.rm = T),
            sd = sd(por_score, na.rm = T),
            n = n()) %>%
  ungroup()

d_all %>% 
  gather(por_scale, por_score, c(pv_score, por_score)) %>%
  mutate(por_scale = recode(por_scale,
                            "pv_score" = "Porosity Vignettes",
                            "por_score" = "Porosity Scale")) %>%
  filter(!is.na(por_score)) %>%
  ggplot(aes(x = country, y = por_score, 
             color = country, fill = country,
             group = religion)) +
  facet_grid(cols = vars(study, por_scale), scales = "free", space = "free") +
  geom_point(position = position_jitterdodge(jitter.width = 0.8,
                                             jitter.height = 0.02,
                                             dodge.width = 0.75), 
             alpha = 0.15) +
  geom_pointrange(data = d_por_sum,
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd, 
                      shape = religion), 
                  position = position_dodge(width = 0.75),
                  fill = "black",
                  color = "black") +
  geom_text(data = d_por_sum %>%
              mutate(ypos = case_when(
                grepl("charismatic", tolower(religion)) ~ mean + sd + 0.05,
                TRUE ~ mean - sd - 0.05)),
            aes(y = ypos, label = paste0("n=", n)), 
            position = position_dodge(width = 0.75),
            size = 3, color = "black") +
  scale_shape_manual(values = 21:24) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  guides(color = F, fill = F, 
         shape = guide_legend(override.aes = list(fill = "black"))) +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(x = "Country", y = "Porosity score (rescaled to 0-1)",
       # caption = "Error bars are ±1 standard deviation from the mean",
       shape = "Religion")

Absorption scores

d_abs_sum <- d_all %>%
  filter(!is.na(abs_score)) %>%
  group_by(study, country, religion) %>%
  summarise(mean = mean(abs_score, na.rm = T),
            sd = sd(abs_score, na.rm = T),
            n = n()) %>%
  ungroup()

d_all %>% 
  filter(!is.na(abs_score)) %>%
  ggplot(aes(x = country, y = abs_score, 
             color = country, fill = country,
             group = religion)) +
  facet_grid(. ~ study, scales = "free", space = "free") +
  geom_point(position = position_jitterdodge(jitter.width = 0.8,
                                             jitter.height = 0.02,
                                             dodge.width = 0.75), 
             alpha = 0.15) +
  geom_pointrange(data = d_abs_sum,
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd, 
                      shape = religion), 
                  position = position_dodge(width = 0.75),
                  fill = "black",
                  color = "black") +
  geom_text(data = d_abs_sum %>%
              mutate(ypos = case_when(
                grepl("charismatic", tolower(religion)) ~ mean + sd + 0.05,
                TRUE ~ mean - sd - 0.05)),
            aes(y = ypos, label = paste0("n=", n)), 
            position = position_dodge(width = 0.75),
            size = 3, color = "black") +
  scale_shape_manual(values = 21:24) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  guides(color = F, fill = F, 
         shape = guide_legend(override.aes = list(fill = "black"))) +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(x = "Country", y = "Absorption score (rescaled to 0-1)",
       # caption = "Error bars are ±1 standard deviation from the mean",
       shape = "Religion")

r1_spev <- lm(spirit_score_std ~ country, d1)
r1_pv <- lm(por_score_std ~ country, d1)
r1_abs <- lm(abs_score_std ~ country, d1)

r2_spev <- lm(spev_score_std ~ country, d2)
r2_dse <- lm(dse_score_std ~ country, d2)
r2_abs <- lm(abs_score_std ~ country, d2)

r3_spev <- lm(spirit_score_std ~ epi_ctry, d3)
r3_por <- lm(por_score_std ~ epi_ctry, d3)

r4_spev <- lm(spev_score_std ~ p7_ctry, d4)
r4_dse <- lm(dse_score_std ~ p7_ctry, d4)
r4_dse <- lm(dse_score_std ~ p7_ctry, d4)
r4_por <- lm(por_score_std ~ p7_ctry, d4)
r4_pv <- lm(pv_score_std ~ p7_ctry, d4)
r4_abs <- lm(abs_score_std ~ p7_ctry, d4)

var	scale	study 1	study 2	study 3	study 4
spiritual experience	DSE	.	.	50%	50%
	spiritual events	22%	24%	38%	31%
porosity	porosity scale	.	61%	.	39%
	porosity vignettes	37%	.	.	27%
absorption	absorption	10%	.	11%	9%

r1_spev_pv <- lm(spirit_score_std ~ por_score_std, d1)
r1_spev_abs <- lm(spirit_score_std ~ abs_score_std, d1)

r2_spev_abs <- lm(spev_score_std ~ abs_score_std, d2)
r2_dse_abs <- lm(dse_score_std ~ abs_score_std, d2)

r3_spev_por <- lm(spirit_score_std ~ por_score_std, d3)

r4_spev_por <- lm(spev_score_std ~ por_score_std, d4)
r4_dse_por <- lm(dse_score_std ~ por_score_std, d4)
r4_spev_pv <- lm(spev_score_std ~ pv_score_std, d4)
r4_dse_pv <- lm(dse_score_std ~ pv_score_std, d4)
r4_spev_abs <- lm(spev_score_std ~ abs_score_std, d4)
r4_dse_abs <- lm(dse_score_std ~ abs_score_std, d4)

df_allb %>% 
  full_join(df_all %>% 
              filter(var == "spiritual experience") %>%
              select(-var) %>%
              rename(outcome = scale) %>%
              mutate(predictor = "country")) %>%
  mutate(predictor = factor(predictor,
                            levels = c("country", "porosity vignettes",
                                       "porosity scale", "absorption"))) %>%
  mutate(percent_exp = paste0(round(rsq * 100), "%")) %>%
  select(-rsq) %>%
  spread(study, percent_exp) %>%
  mutate_at(vars(starts_with("study")),
            funs(case_when(is.na(.) ~ ".", 
                           TRUE ~ .))) %>%
  kable() %>% 
  kable_styling() %>%
  collapse_rows(1:2)

Joining, by = c("predictor", "outcome", "study", "rsq")
Column `predictor` joining factor and character vector, coercing into character vectorColumn `outcome` joining factor and character vector, coercing into character vector

predictor	outcome	study 1	study 2	study 3	study 4
country	DSE	.	.	50%	50%
	spiritual events	22%	24%	38%	31%
porosity vignettes	DSE	.	.	.	28%
	spiritual events	13%	.	.	29%
porosity scale	DSE	.	.	.	44%
	spiritual events	.	39%	.	33%
absorption	DSE	.	.	3%	2%
	spiritual events	11%	.	12%	8%

Relationships

All studies, multipart plot

Spiritual Events only, cowplot

fig_s1_por <- d_all %>%
  filter(study == "Study 1") %>%
  ggplot(aes(x = pv_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Vignettes",
       y = "Spiritual Events",
       color = "Country")
# fig_s1_por

fig_s1_abs <- d_all %>%
  filter(study == "Study 1") %>%
  ggplot(aes(x = abs_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Absorption",
       y = "Spiritual Events",
       color = "Country")
# fig_s1_abs

fig_s1_title <- ggdraw() + 
  draw_label("STUDY 1", fontface = 'bold', x = 0, hjust = 0.5) +
  theme(plot.margin = margin(0, 0, 0, 105))

fig_s1 <- plot_grid(
  fig_s1_title,
  plot_grid(fig_s1_por, fig_s1_abs, ncol = 1, labels = c("A", "B")),
  ncol = 1, rel_heights = c(1, 10))

Removed 20 rows containing non-finite values (stat_smooth).Removed 20 rows containing non-finite values (stat_smooth).Removed 20 rows containing missing values (geom_point).Removed 28 rows containing non-finite values (stat_smooth).Removed 28 rows containing non-finite values (stat_smooth).Removed 28 rows containing missing values (geom_point).

# fig_s1

fig_s1_title_vert <- ggdraw() + 
  draw_label("STUDY 1", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s1_vert <- plot_grid(
  fig_s1_title_vert,
  plot_grid(fig_s1_por, fig_s1_abs, ncol = 2, labels = c("A", "B")),
  ncol = 1, rel_heights = c(1, 10))

Removed 20 rows containing non-finite values (stat_smooth).Removed 20 rows containing non-finite values (stat_smooth).Removed 20 rows containing missing values (geom_point).Removed 28 rows containing non-finite values (stat_smooth).Removed 28 rows containing non-finite values (stat_smooth).Removed 28 rows containing missing values (geom_point).

# fig_s1_vert

fig_s2_abs <- d_all %>%
  filter(study == "Study 2") %>%
  ggplot(aes(x = abs_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Absorption",
       y = "Spiritual Events",
       color = "Country")
# fig_s2_abs

fig_s2_title <- ggdraw() + 
  draw_label("STUDY 3", fontface = 'bold', x = 0, hjust = 0.5) +
  theme(plot.margin = margin(0, 0, 0, 105))

fig_s2 <- plot_grid(
  fig_s2_title,
  plot_grid(NULL, fig_s2_abs, ncol = 1, labels = c("", "C")),
  ncol = 1, rel_heights = c(1, 10))

Removed 10 rows containing missing values (geom_smooth).

# fig_s2

fig_s2_title_vert <- ggdraw() + 
  draw_label("STUDY 3", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s2_vert <- plot_grid(
  fig_s2_title_vert,
  plot_grid(NULL, fig_s2_abs, ncol = 2, labels = c("", "C")),
  ncol = 1, rel_heights = c(1, 10))

Removed 10 rows containing missing values (geom_smooth).

# fig_s2_vert

fig_s3_por <- d_all %>%
  filter(study == "Study 3") %>%
  ggplot(aes(x = por_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Scale",
       y = "Spiritual Events",
       color = "Country")
# fig_s3_por

fig_s3_title <- ggdraw() + 
  draw_label("STUDY 2", fontface = 'bold', x = 0, hjust = 0.5) +
  theme(plot.margin = margin(0, 0, 0, 105))

fig_s3 <- plot_grid(
  fig_s3_title,
  plot_grid(fig_s3_por, NULL, ncol = 1, labels = c("D", "")),
  ncol = 1, rel_heights = c(1, 10))
# fig_s3

fig_s3_title_vert <- ggdraw() + 
  draw_label("STUDY 2", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s3_vert <- plot_grid(
  fig_s3_title_vert,
  plot_grid(fig_s3_por, NULL, ncol = 2, labels = c("D", "")),
  ncol = 1, rel_heights = c(1, 10))
# fig_s3_vert

fig_s32_vert <- plot_grid(
  plot_grid(
    fig_s3_title_vert, 
    plot_grid(fig_s3_por, labels = c("C")), 
    ncol = 1, rel_heights = c(1, 10)),
  plot_grid(
    fig_s2_title_vert,
    plot_grid(fig_s2_abs, labels = c("D")), 
    ncol = 1, rel_heights = c(1, 10))
)

Removed 10 rows containing missing values (geom_smooth).

fig_s4_por1 <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = pv_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Vignettes",
       y = "Spiritual Events",
       color = "Country")
# fig_s4_por1

fig_s4_por2 <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = por_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Scale",
       y = "Spiritual Events",
       color = "Country")
# fig_s4_por2

fig_s4_abs <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = abs_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Absorption",
       y = "Spiritual Events",
       color = "Country")
# fig_s4_abs

fig_s4_title <- ggdraw() + 
  draw_label("STUDY 4", fontface = 'bold', x = 0, hjust = 0.5) +
  theme(plot.margin = margin(0, 0, 0, 190))

fig_s4 <- plot_grid(
  fig_s4_title,
  plot_grid(plot_grid(fig_s4_por1, fig_s4_por2, ncol = 2, labels = c("E", "F")), 
            plot_grid(NULL, fig_s4_abs, NULL, ncol = 3, rel_widths = c(1, 2, 1), labels = c("", "G", "")), 
            ncol = 1),
  ncol = 1, rel_heights = c(1, 10))

Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).Removed 11 rows containing missing values (geom_smooth).

# fig_s4

fig_s4_title_vert <- ggdraw() + 
  draw_label("STUDY 4", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s4_vert <- plot_grid(
  fig_s4_title_vert,
  plot_grid(plot_grid(fig_s4_por1, fig_s4_por2, ncol = 1, labels = c("E", "F")), 
            plot_grid(NULL, fig_s4_abs, NULL, ncol = 1, rel_heights = c(1, 2, 1), labels = c("", "G", "")), 
            ncol = 2),
  ncol = 1, rel_heights = c(1, 20))

Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).Removed 11 rows containing missing values (geom_smooth).

# fig_s4_vert

fig_legend <- get_legend(fig_s4_por1 + theme(legend.position = "bottom"))

Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).

fig_all <- plot_grid(fig_s1, fig_s2, fig_s3, fig_s4, ncol = 4,
                     rel_widths = c(1, 1, 1, 2), scale = 0.95)
fig_all

fig_all_vert <- plot_grid(fig_s1_vert, fig_s32_vert, fig_s4_vert, fig_legend,
                          ncol = 1, rel_heights = c(1, 1, 2, 0.2))
fig_all_vert

Daily Spiritual Experiences only, cowplot

fig_s2_abs <- d_all %>%
  filter(study == "Study 2") %>%
  ggplot(aes(x = abs_score, y = dse_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Absorption",
       y = "Daily Spiritual Experiences",
       color = "Country")
# fig_s2_abs

fig_s2_title_vert <- ggdraw() + 
  draw_label("STUDY 3", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s2_vert <- plot_grid(
  fig_s2_title_vert,
  plot_grid(fig_s2_abs, ncol = 1, labels = "B"),
  ncol = 1, rel_heights = c(1, 10))
# fig_s2_vert

fig_s4_por1 <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = pv_score, y = dse_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Vignettes",
       y = "Daily Spiritual Experiences",
       color = "Country")
# fig_s4_por1

fig_s4_por2 <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = por_score, y = dse_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Scale",
       y = "Daily Spiritual Experiences",
       color = "Country")
# fig_s4_por2

fig_s4_abs <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = abs_score, y = dse_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Absorption",
       y = "Daily Spiritual Experiences",
       color = "Country")
# fig_s4_abs

fig_s4_title_vert <- ggdraw() + 
  draw_label("STUDY 4", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s4_por1_vert <- plot_grid(
  fig_s4_title_vert,
  plot_grid(fig_s4_por1, ncol = 1, labels = "A"),
  ncol = 1, rel_heights = c(1, 10))
# fig_s4_por1_vert

fig_s4_por2_vert <- plot_grid(
  fig_s4_title_vert,
  plot_grid(fig_s4_por2, ncol = 1, labels = "C"),
  ncol = 1, rel_heights = c(1, 10))
# fig_s4_por2_vert

fig_s4_abs_vert <- plot_grid(
  fig_s4_title_vert,
  plot_grid(fig_s4_abs, ncol = 1, labels = "D"),
  ncol = 1, rel_heights = c(1, 10))
# fig_s4_abs_vert

fig_legend <- get_legend(fig_s4_por1 + theme(legend.position = "bottom"))

fig_all_vert <- plot_grid(
  plot_grid(fig_s4_por1_vert, fig_s2_vert, 
            fig_s4_por2_vert, fig_s4_abs_vert,
            ncol = 2),
  fig_legend,
  ncol = 1, rel_heights = c(2, 0.2))
fig_all_vert

Other versions

Spiritual Events only, one grid, new layout

d_all %>%
  gather(spirit_scale, spirit_score, c(spev_score, dse_score)) %>%
  gather(pred_scale, pred_score, c(por_score, pv_score, abs_score)) %>%
  mutate(spirit_scale = recode_factor(spirit_scale,
                                      "spev_score" = "Spiritual Events",
                                      "dse_score" = "Daily Spiritual Experiences"),
         pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption"),
         pred_type = case_when(pred_scale == "Absorption" ~ "Absorption",
                               grepl("Porosity", pred_scale) ~ "Porosity",
                               TRUE ~ NA_character_),
         pred_type = factor(pred_type, levels = c("Porosity", "Absorption")),
         study_scale = paste(study, pred_scale, sep = ": "),
         study_scale2 = case_when(
           study == "Study 4"  & pred_scale != "Absorption" ~ 
             paste(study, pred_scale, sep = ": "),
           TRUE ~ study),
         study_scale3 = case_when(
           study == "Study 4" & pred_scale == "Porosity Scale" ~ "Porosity Scale",
           study %in% c("Study 1", "Study 4") ~ "Porosity Vignettes",
           study == "Study 2" ~ "Porosity Scale", 
           TRUE ~ " "),
         study_scale3 = factor(study_scale3,
                               levels = c("Porosity Vignettes",
                                          "Porosity Scale", " "))) %>%
  filter(!is.na(pred_score),
         spirit_scale == "Spiritual Events") %>%
  ggplot(aes(x = pred_score, y = spirit_score)) +
  facet_grid(rows = vars(pred_type), cols = vars(study, study_scale3)) +
  # facet_grid(pred_type ~ study_scale3) +
  geom_point(data = . %>% distinct(study, study_scale3, country, 
                                   pred_type, pred_scale, pred_score,
                                   spirit_scale, spirit_score),
             aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Score on predictor scale (Porosity Vignettes, Porosity Scale, or Absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

Spiritual Events only, by study

fig_s1 <- d_all %>%
  filter(study == "Study 1") %>%
  distinct() %>%
  gather(pred_scale, pred_score, c(pv_score, por_score, abs_score)) %>%
  mutate(pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption")) %>%
  # mutate(study_scale = paste(study, pred_scale, sep = ": ")) %>%
  filter(!is.na(pred_score)) %>%
  ggplot(aes(x = pred_score, y = spev_score)) +
  facet_grid(pred_scale ~ study) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm",
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  # theme(legend.position = "bottom", 
  #       axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  # guides(color = guide_legend(override.aes = list(alpha = 1))) +
  guides(color = F) +
  labs(x = "Score on predictor scale\n(porosity or absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events scale (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

fig_s2 <- d_all %>%
  filter(study == "Study 2") %>%
  distinct() %>%
  gather(pred_scale, pred_score, c(pv_score, por_score, abs_score)) %>%
  mutate(pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption")) %>%
  # mutate(study_scale = paste(study, pred_scale, sep = ": ")) %>%
  filter(!is.na(pred_score)) %>%
  ggplot(aes(x = pred_score, y = spev_score)) +
  facet_grid(pred_scale ~ study) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm",
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  # theme(legend.position = "bottom", 
  #       axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  # guides(color = guide_legend(override.aes = list(alpha = 1))) +
  guides(color = F) +
  labs(x = "Score on predictor scale\n(porosity or absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events scale (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

fig_s3 <- d_all %>%
  filter(study == "Study 3") %>%
  distinct() %>%
  gather(pred_scale, pred_score, c(pv_score, por_score, abs_score)) %>%
  mutate(pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption")) %>%
  # mutate(study_scale = paste(study, pred_scale, sep = ": ")) %>%
  filter(!is.na(pred_score)) %>%
  ggplot(aes(x = pred_score, y = spev_score)) +
  facet_grid(pred_scale ~ study) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm",
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  # theme(legend.position = "bottom", 
  #       axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  # guides(color = guide_legend(override.aes = list(alpha = 1))) +
  guides(color = F) +
  labs(x = "Score on predictor scale\n(porosity or absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events scale (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

fig_s4 <- d_all %>%
  filter(study == "Study 4") %>%
  distinct() %>%
  gather(pred_scale, pred_score, c(pv_score, por_score, abs_score)) %>%
  mutate(pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption")) %>%
  # mutate(study_scale = paste(study, pred_scale, sep = ": ")) %>%
  filter(!is.na(pred_score)) %>%
  ggplot(aes(x = pred_score, y = spev_score)) +
  facet_grid(pred_scale ~ study) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm",
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  # theme(legend.position = "bottom", 
  #       axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  # guides(color = guide_legend(override.aes = list(alpha = 1))) +
  guides(color = F) +
  labs(x = "Score on predictor scale\n(porosity or absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events scale (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

Spiritual Events only, full grid

d_all %>%
  distinct() %>%
  gather(pred_scale, pred_score, c(pv_score, por_score, abs_score)) %>%
  mutate(pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption")) %>%
  # mutate(study_scale = paste(study, pred_scale, sep = ": ")) %>%
  filter(!is.na(pred_score)) %>%
  ggplot(aes(x = pred_score, y = spev_score)) +
  facet_grid(pred_scale ~ study) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm",
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  # theme(legend.position = "bottom", 
  #       axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  # guides(color = guide_legend(override.aes = list(alpha = 1))) +
  guides(color = F) +
  labs(x = "Score on predictor scale\n(porosity or absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events scale (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

Broken down by predictor type and spiritual scale

d_all %>%
  gather(spirit_scale, spirit_score, c(spev_score, dse_score)) %>%
  gather(poros_scale, poros_score, c(por_score, pv_score)) %>%
  mutate(spirit_scale = recode_factor(spirit_scale,
                                      "spev_score" = "Spiritual Events",
                                      "dse_score" = "Daily Spiritual Experiences"),
         poros_scale = recode_factor(poros_scale,
                                     "pv_score" = "Porosity Vignettes",
                                     "por_score" = "Porosity Scale"),
         study_scale = paste(study, poros_scale, sep = ": ")) %>%
  filter(!is.na(poros_score)) %>%
  ggplot(aes(x = poros_score, y = spirit_score)) +
  facet_grid(spirit_scale ~ study_scale) +
  geom_point(data = . %>% distinct(study, study_scale, country, 
                                   poros_scale, poros_score,
                                   spirit_scale, spirit_score),
             aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Score on porosity measure (rescaled to 0-1)",
       y = "Score on spiritual experience measure (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

d_all %>%
  gather(spirit_scale, spirit_score, c(spev_score, dse_score)) %>%
  mutate(spirit_scale = recode_factor(spirit_scale,
                                      "spev_score" = "Spiritual Events",
                                      "dse_score" = "Daily Spiritual Experiences"),
         study_scale = paste(study, "Absorption scale", sep = ": ")) %>%
  filter(!is.na(abs_score)) %>%
  ggplot(aes(x = abs_score, y = spirit_score)) +
  facet_grid(spirit_scale ~ study_scale) +
  geom_point(data = . %>% distinct(study, study_scale, country, 
                                   abs_score,
                                   spirit_scale, spirit_score),
             aes(color = country), alpha = 0.2) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Score on absorption measure (rescaled to 0-1)",
       y = "Score on spiritual experience measure (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")